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Abstract In this paper we propose methods for computing Fresnel integrals 
based on truncated trapezium rule approximations to integrals on the real 
line, these trapezium rules modified to take into account poles of the integrand 
near the real axis. Our starting point is a method for computation of the error 
function of complex argument due to Matta and Reichel (J. Math. Phys. 34 
(1956), 298-307) and Hunter and Regan (Math. Comp. 26 (1972), 539-541). 
pH We construct approximations which we prove are exponentially convergent as 

a function of N, the number of quadrature points, obtaining explicit error 
bounds which show that accuracies of 10~ 15 uniformly on the real line are 
achieved with N = 12, this confirmed by computations. The approximations 
we obtain are attractive, additionally, in that they maintain small relative 
errors for small and large argument, are analytic on the real axis (echoing the 
analyticity of the Fresnel integrals), and are straightforward to implement. In 
a last section we explore the implications of our results for the computation 
, of the error function of real and complex argument. 
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1 Introduction 

Let C(x), S(x), and F(x) be the Fresnel integrals denned by 

C(x):= cos (±7rf 2 ) dt, S(x):= sin (^i 2 ) di, (1) 
Jo J a 

and 

p-wr/4 /-oo i 

F(s) := ?—=- / e" dt. (2) 



Our definitions in ([T]) are those of [2] and [TJ §7.2(iii)], and F, C and S are 
related through 

y/2e™' A F(x) = \-C (y/2/irx\ -5 (V^r*)) ■ ( 3 ) 

Another definition in common use [TJ §7.2(iii)] is 

/•oo 

T(x):=\ e^ it2 dt = V2c i7T / 4 F(^j2x). (4) 

Fresnel integrals arise in applications throughout science and engineering, 
especially in problems of wave diffraction and scattering (e.g., [H §8.2], [5]), so 
that methods for the efficient and accurate computation of these functions are 
of wide application. The purpose of this paper is to present new approximations 
for the Fresnel integrals, based on iV-point trapezium rule approximations to 
integral representations for F(x), these trapezium rules modified to take into 
account the presence of poles in the integrand near the path of integration. 
The resulting approximation to F(x) is 

*<x) := i + I^W") £^ (5) 

k— 1 

N 



X 



3 i(x^+7r/4) 



E^2-^2. ( 6 ) 



where 



exp (2A N xe- i7T / 4 ) + 1 A N k~i x2+it l 



(k - 1/2) 7T 



tk:= )==y = > A w :=ijv+i = V(^ + 1/2)tt. (7) 

v (,-<* + V 2 )^ 

The corresponding approximations to C(x) and S(x) (obtained by substituting 
this approximation in ([3]) and separating real and imaginary parts) are 

£, / n 1 sinh (t/tt An x) +sin(y / 7rAjva;) 

2 cos( v / 7r Ajv x) + cosh( v / 7r^47v x) 

a N (|a; 2 ) sin (|a; 2 ) - & w (|x 2 ) cos (|a; 2 ) ) (8) 



.4 
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and 

S N (x) 
where 



1 sinh (y/w An x) — sin (y/n An x) 

2 cos(-v/7r An x) + cosh(y / 7r Am x) 

a N (s) :=*E^T74' M«):=£^^T- ( 10 ) 

fc=i k fe=i fc 

These approximations, designed for computation of F(x), C{x) and S(x) 
for all x G M, are attractive in several respects. Firstly, they are provably 
exponentially convergent as N increases: a main result of the paper is to show 
in ^2] that, for N G N, where 

E N (x) :=F{x) -F N {x), (11) 

we have 

— ttN — 7T N 

\E N (x)\ < c N ; e < c (12) 

y/N+1/2 y/N + 1/2 

Here cjv is a decreasing sequence of positive constants, c\ > c-2 > given 
explicitly by 

20V2e-' r / 2 / „ _ flA »\ (2tt + lje"^ 2 

C N = -T- 1 + 2y/nc p An + - — - — '- , (13) 

97r(l-e- 2A «) V V J 2y/2^/ 2 A N V ; 

where (3 ~ 0.0536 is given by (|4"2")l . so that 

20v / 2e~ ,r/2 

d w 0.825 and lim C7V = — « 0.208. (14) 

N-S-oo 97T 



The bound (|12|) implies exponential convergence also of the approximations 
for C(x) and S(x). Indeed, since holds with F, C, and S replaced by Fn, 
Cn, and Sn, respectively, it holds for x G K (sec SJ3]) that 

\C(x)-C N (x)\ < V2\E N {y/^x)\ and \S(x)-S N (x)\ < \/2\E N (y/^/2x)\. 

(15) 

We will present in 21 numerical computations which demonstrate this ex- 
ponential convergence and show that, with only N = 12 quadrature points, 
the absolute error \En(x)\ < 10" 15 and the same bound holds for the relative 
error, i.e., \En(x)\ < 10 _15 |F(x)|, both these bounds holding for all x G R. 
Thus these approximations are efficient; they achieve close to double precision 
machine accuracy with very little computation. 

A second attractive feature of these approximations is their smoothness. 
The Fresnel integrals F(z), C(z), and S(z), considered as functions of z G C in 
the complex plane, are all entire functions. The approximations Fn{z), Cn{z) 
and Sn{z) are meromorphic rather than entire functions, but are analytic in 
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a strip which surrounds the real axis, in particular arc infinitely diffcrcntiable 
on the real line. But more than this, the bound (|T2"]) extends into the com- 
plex plane, holding in the first and third quadrants, and, in modified form 
(see ([57]) ), in the strip |Im(z)| < An/(2^/2) around the real axis. This implies 
exponentially convergent error estimates, presented in Sj2]and for the differ- 
ence between the coefficients in the Maclaurin series of F, C, and S and those 
in the corresponding series for Fjy, Cm and Sn- In turn (see this implies 
that the approximations are particularly accurate and retain small relative 
error for \x\ small, and the computations in 2] demonstrate this. 

The third attractive feature of these approximations is that they inherit 
certain symmetries of the Fresnel integrals. In particular, our normalisation of 
F(x) is such that 

F(-x) = 1 - F(x), (16) 

so that, in particular, F(0) = 1/2. It is clear from ([5]) that the same holds for 
F N (x), i.e., 

F N (-x) = l-F N (x). (17) 
Similarly, where an overline denotes a complex conjugate, 

F(z) = F(iz) and F N (z) = F N (iz). (18) 

Both these symmetries can be deduced as a consequence of the structure of C 
and S and their approximations: by inspection of ([8]) and ([9]) we see that 

C N {x) = xf c (x 4 ), S N (x) = x 3 f s (x 4 ), (19) 

where fc and fs are analytic in a neighbourhood of the real line and are real- 
valued for real arguments. This is the same structure as C and S (see |67|) 
below). In particular, (fT9"|) implies that Cm and Sn, like C and S, are odd 
functions. 

The final attractive feature is that these approximations are straightfor- 
ward to code. Tabic Q] shows the simple Matlab code used to evaluate Fm(x) 
for all the computations in this paper. Of course this code is easily converted 
to other computer languages. 

Let us make some comments regarding the antecedents of our method. 
The methodology we employ has a long history, though we introduce in this 
paper significant improvements in implementation and in error analysis. The 
derivation of our approximation Fm (x) makes use of the relationship between 
the Fresnel integral and the error function, that 

F(x) = icrfc(c~ i7r/4 x) = | c ix2 w (c i7r / 4 a;) (20) 

where crfc is the complementary error function, defined by 



erfc(z) := 




and 

2 

w(z) := e~ z erfc(— iz). 
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function f = fresnel(x,N) 

"/, Evaluates the approximation F_N(x) to the Fresnel integral F(x). 
% x is a real scalar or matrix, 

% N is the positive integer controlling accuracy (suggest N=12) , 
% f is the corresponding scalar or matrix of values of F_N(x) . 
f = zeros (size (x) ) ; 
select = x>=0; 

if any (select), f (select) = F(x(select) ,N) ; end 

if any(~select) , f ("select) = l-F(-x(~select) ,N) ; end 

function f = F(x,N) 

h = sqrt(pi/(N+0.5)) ; 

t = h*((N:-l:l)-0.5); AN = pi/h; 

t2 = t.*t; t4 = t2.*t2; et2 = exp(-t2) ; 

rooti = exp(i*pi/4); 

z = rooti*x; x2 = x.*x; x4 = x2.*x2; z2 = i*x2; 
S = (-et2(l) ./(x4+t4(l))) .*(z2+t2(l)) ; 
for n = 2:N 

S = S + (-et2(n) . / (x4+t4(n) ) ) .*(z2+t2(n)) ; 

end 

ez = exp((2*AN*i*rooti)*x) ; 

f = (i/AN)*z.*exp(z2) . *S + ez./(ez+l); 



Table 1 Matlab code to evaluate Fjsr(x) given by JJj}, making use of 11171 for x < 0. 



It also depends on the integral representation J5J (7.1.4)] that 

i r°° e~* 2 iz f°° e"* 2 

w(z) = - / dt=- — - dt, for Im(z) > 0. (21) 

7T J_ OQ z-t 7T J_ ao z 2 - t- 

Combining (|2Tj|) and ([2"Tj) gives an integral representation for F(x), that 

F(x) = — S x2+7T/4) ( 6 * „ dt, for x > 0. (22) 

The observation that the trapezium rule is exponentially convergent when 
applied to integrals of the form 

/oo 
f(t)e~ t2 dt, (23) 

with f(t) analytic in a strip surrounding the real axis, dates back at least to 
Turing |28j and Goodwin |12j . The derivation of this result uses contour inte- 
gration and Cauchy's residue theorem; see $5] below. Applying the trapezium 
rule with step-length ft, > to (f2"2"|) leads to the approximation 

OO _ T '2 , OO _ r 2 

k= — oo K k=l K 

(24) 

where 

r fc := (k - l/2)h. (25) 
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When x > is large this approximation is very accurate, indeed is essentially 
identical to the approximation Fn (x) with ./V large if the choice 



h = tt/A n = y/n/(N + l/2) (26) 

is made. However, this approximation becomes increasingly poor as x > 
approaches zero. 

In the context of developing methods for evaluating the complementary 
error function of complex argument (by ([27J]) . evaluating F(x) for x real is just 
a special case of this larger problem), Chiarella and Riechel [7], Matta and 
Reichel [18] , and Hunter and Regan [15] proposed modifications of the trapez- 
ium rule that follow naturally from the contour integration argument used to 
prove that the trapezium rule is exponentially convergent. The most appro- 
priate form of this modification is that in j!5j where the modified trapezium 
rule approximation 

F(x) pa — e i <" 2 +"/ 4 ) 2 ■ 2 +R(h,x), for x > 0, (27) 
n k=i X 1Tfc 

is proposed. Here the correction term R(h,x) is defined by 

C l/(exp(27rc- i,r /V^) + if < x < y/^ir/h, 
R{h,x):= < 0.5/(exp(27rc- i7r / 4 a;//i) + 1), if x = V2ir/h, 
[o, tfx>V2n/h. 

The approximation ((27)) clearly coincides with Fjv(x) for < x < \f2-Kjh if 
the range of summation in ([2T]) is truncated to 1, N and the choice (|2"?H for 
h is made. Hunter and Regan prove that the magnitude of the error in the 
approximation (1271) is 

xe -ir 2 /h 2 

~ 0F(l-e~ 2 - 2 A 2 ) |x 2 /2-7r 2 //i 2 |' 

for x > 0, provided x ^ \f2i:/h. Similar estimates, it appears arrived at 
independently are derived by Mori [19] . in which paper the emphasis is on 
computing erfc(x) for real x. 

It may be becoming clearer to the reader at this point what the contribu- 
tions of this paper are, which are to improve and develop the work of Hunter 
and Regan |15| and the articles that precede it, as they apply to the compu- 
tation of Fresnel integrals. The new contributions in this paper are essentially 
four. The first is to truncate the range of summation in ([2"Tf and to propose the 
choice (|26p in order to balance the error arising from truncating the infinite 
sum to the finite range 1, ...,N with the error made in the trapezium rule ap- 
proximation modified in the spirit of |15| . With this change the approximation 
([2"T]) becomes identical to Fn(x) for < x < ^/2-K/h. The second contribu- 
tion is to recognise that (|2"T|) . while accurate, has the drawback that R(x, h) is 
discontinuous, so that the entire function F(x) is being approximated on the 
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real line by a discontinuous function (albeit with small discontinuities). This 
contribution is further to realise that the approximation formula proposed on 
< x < ^/2n/h in fact provides a smooth and accurate approximation on the 
whole real line. The third contribution, the contribution which is most sub- 
stantial in terms of analysis, is to improve the error bound (|28p of [T5]. This 
error bound is unsatisfactory in that it blows up at x = \f2i{/h in a way not 
seen in the numerical results in SJH (though there is some increase observed in 
the error when x is near this value; see Figure [3]). The bounds we prove in Sj4] 
modify the arguments of [IB] and [TS] to establish rigorous bounds, on both 
the absolute and relative error in the approximation Fn(x), that are uniform 
in x on the whole real line. (These rigorous bounds include, in particular, the 
bound The final contribution is to explore the implications of these re- 

sults for the computation of C(x) and S(x), and, briefly, for the computation 
of erfc(x) for real and complex argument x. 

As mentioned above, it is clear from ([20]) that numerical methods devel- 
oped to evaluate erfc(z) for complex z can be applied in particular to compute 
F(x) for x £ R. Indeed the starting point for the new method we propose 
is the method in [TS] for erfc(z). Since [T5] was published, further algorithms 
for computing erfc(z) with complex argument z have been developed. The 
algorithm having the best combination of efficiency and accuracy for inter- 
mediate values of z (say in the range 1.5 < \z\ < 5), where computation is 
most challenging, appears to be that of Weideman [55]. We will make compar- 
ison with this algorithm in f[¥] showing that our algorithm appears to achieve 
substantially higher accuracy with fewer computations. 

Of course, given the importance of Fresnel integrals in applications, there 
exist also many effective methods specifically designed to evaluate C(x), S(x) 
and/or F(x) for real x. We will describe the best of these in ^3.11 though 
we note that none of these approximations appear to have the combination 
of features of our new approximations ([5]), ©, and ©, namely that they 
are analytic in a strip around the real line, are exponentially convergent with 
provable error bounds, and are straightforward to implement. 

We end this introduction by outlining the remainder of the paper. In 3Hwc 
derive the approximation ([5]) to F(x) and prove rigorous bounds on En(x), 
including ([T2"]) . In S[3]we deduce from this the approximations and © and 
bounds on the errors C(x) — Cn(x) and S(x) — Sn{x), especially bounds for x 
small, and survey other methods for computing Fresnel integrals. In 2)we show 
numerical results, comparing our new approximations with the error bounds 
derived in the earlier sections and with certain rival methods for computing 
Fresnel integrals. In concluding remarks in iJS]we explain how the methods we 
introduce in Sj2] are potentially of much wider application to computation of 
special functions, not least erfc(z) for complex z and other functions arising in 
computational acoustics. The appendix proves a sharp lower bound for |-F(x)| 
for x > 0, of independent interest and a key component in our theoretical 
bounds on relative errors. Maintaining a theme, the proof of this lower bound 
requires modified trazezium rule approximations and their error estimates. 
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2 The approximation of F(x) and its error bounds 

In this section we derive the approximation ([3]) to F(x) and derive (|12l) and 
related error bounds for this approximation, these error bounds demonstrating 
that both the absolute and relative errors in the approximation Fn (x) converge 
exponentially to zero as N increases, uniformly on the real line, and that 
N = 12 is large enough to achieve errors < 10 -15 . 

The first part of our derivation follows in large part Matta and Reichel [T8] 
and Hunter and Regan [TS]- From (|2"2"|) we have that, for a; > 0, 

J := / /(*) dt = F(x), where /(f) := e^ f - 5 , (29) 

J_ OQ 2tt x z + it z 

and we have suppressed in our notation the dependence of /(f) on x. Choose 
a step-size h > for the trapezium rule and let 

g{z) = itan(7rz//i), 

which is a meromorphic function with simple poles at the points r^, denned 
by ([25]), which has the property that, for z = X + \H with X e E, H > 0, 

|l+g(z)[< 1 _ e _ 27rg/fe . (30) 

The approximation (|27[) is obtained by considering the integral in the complex 
plane, 

J = J f(z){l + g(z))dz, 

where the path of integration is from — oo to oo along the real axis, except 
that the path makes small semicircular deformations to pass above each of the 
simple poles at the points Tfc, k £ 1. Explicitly, the fcth deformation is the 
semicircle jk = {tr + ee _l61 : n < 6 < 2ir}, with e in the range (0,/i/2) small 
enough so that the simple pole singularity in f(z) at z = zq := e™/ 4 x lies 
above r. Then, since f(z)g(z) is an odd function, we see that 

J = f f(z)dz+ f f(z)g(z)dz = I + Y / I f{z)9{z)dz. 

In the limit e -> 0, f(z)g(z) dz -> -in Res(/g, r fc ) = -tiffa), where 
Res(g, Tfc) denotes the residue of /g at r^. Thus J = I — Ih, where 

oo 

h = hJ2 f(r k ) = 2h J2 ~ 1/2)^) (31) 

fcGZ k=l 

is a trapezium/midpoint rule approximation to /. On the other hand, where 
Th = {x + iH : x G M}, by the residue theorem, 

J = y /(z)(l+. 9 (z))dz + H(V2i/-a;) PC ft , 
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for H > 0, where H is the Heaviside step function and 

1 + itan (c i7r/4 :E7r//i)) . 



PC h = 27riRcs(/(l + g), zo) = \0- + ff(*>)) = \ 



Thus 



I = I h + H (V2 H - sc) PC h + jf f(z)(l + g{z)) dz 



(32) 



The point here is that the integral over Th can be negligible so that a 
good approximation is obtained by the modified trapezium rule approximation, 
//, + H (V2H — x) PCh- In particular, noting (j30|) and that, for z = X + iH, 

|.x 2 +iz 2 | = \z ~z\\z + z\ > \x/V2-H\\x/V2 + H\ = \x 2 /2-H 2 \ 



and that e * dt = ^/t:, we see that 
f(z)(l + g(z))dz 



L 



< 



X c 



H 2 -2-xH/h 



^\H 2 -x 2 /2\ (l-e- 2 ^/' 1 ) 



Choosing H = ir/h, to minimise the exponent H 2 — 2nH/h, it follows that 
I = I h + H (V2n/h - x) PC h + e h with 



M < 5i(x) 



x e 



^\n 2 /h 2 -x 2 /2\ (1-c- 2 - 2 /^ 2 ) 



(33) 



Note that Ih + H (-^2-k/Ii — x) PCh = hi + R{h,x) is precisely the approxi- 
mation (|27[) . and that the above bound on eh is precisely the bound (|28[) from 

na. 

implies that 



Let I* h :=I h + PC h and e* h := I- !£. Then 



Since, applying (f3"U)) , 



3hl < 6i(x), for < x < V2n/h. 

e -V2TTx/h 



(34) 



\PC h \ < 



we sec that 



\e* h \ < S 3 (x) := ^(z) + 



1 



2— y/2-nx/h 



, for x > V2 7r/h. 



(35) 



The bounds (|3"4"1) and (|3"B"|) both blow up as x approaches \p2ixjh. Contin- 
uing to choose H = ir/h, select e in the range (0,H) and consider the case 
that 

x 



V2 



H 



< e. 



(36) 



In this case we observe that the derivation of (f3"2")l can be modified to show 
that 



f(z)(l + g(z))dz 



(37) 



rs 
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where the contour passes above the pole in / at z$; precisely, is the 
union of r' and 7, where r' = {z e Th ■ \z — Zo| > e } an d 7 is the circular 
arc 7 = {z + ec ie : 9 < 9 < ir - 9 }, where 9 = sin _1 ((# - x/Vty/e) G 
(-- 7r/2,7r/2). For z G T' it holds that 



|x- 2 + iz 2 | = |z„ - z| \zo + z\>e \x/V2 + H\ 
Thus, and applying (j3T))) . similarly to (fM]) we deduce that 



/(*)(! + <?(*)) dz 



< 



/TTelTr/Zi + x/v^l (1 -e- 2 * 2 /' 12 ) 



(38) 



(39) 



To bound the integral over 7 we note that, for z = X + \Y = zq + ec iS E 7, 
([3"B")) is true and Y > H. Further, \c~ z \ = e p , where 

P = Y 2 -X 2 = 2xesin(6»-7r/4)-e 2 cos(26l) < 2xe + e 2 < 2\/2i7e + (2\/2 + l)e 2 , 

using (jSHJ). From these bounds and (|3T))) . defining a = e/i? G (0, 1), we deduce 
that 



f(z)(l + g(z))dz 



< 2z exp((2 v / 2a + (2^2 + l)a 2 - 2)7r 2 /fr 2 ) 



s|7r/fc + a?/>/2| (l-e- 2 ^ 2 /^ 2 ) 



Choosing a = 1/4, we can bound e* h using (|3"T|) . (f5U|) . (|4"D")) , and the triangle 
inequality, to get that 



141 < S 2 (x) 



Ahx e'^^ 2 



ir s / 2 \ir/h + x/V2\ (l-e" 2 ^/^) 
for x in the range (|36[) . where 

\/2 2\/2 + 1 



(l + 2V 7 7e-^ 2 /' 12 ) , (41) 



/3 = 1 - 2v / 2a - (2V2 + l)a 2 = 1 - 



16 



0.0536. 



(42) 



Summarising up to this point, we see that we have shown that I = F(x) 
can be approximated by the modified trapezium rule IT = Ih + PCh with error 
e* h = I — 1^. This error satisfies \e* h \ < Ah(x), for x > 0, where 



A h (x) 



Sx(x), 0<^<fff, 
S 2 (x), \H<^< 



(43) 



dSl) , and ([55|1. respec- 



with = 7r//i. Here 5%, 62, and ^3 are defined by 
tively. 

The approximation Fn(x), given by ([5]), that we propose for / = F(x) is 
just Ifa = Ih + PCh with a particular choice of h and with the range of sum- 
mation in pip reduced to the finite range 1, ...,N. This induces an additional 
error, 



T N :=2h J2 /( T « 



(44) 



m=N+l 
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which, for x > 0, is bounded by 

hx e _Tm 

/ 4 T T 4 
-N+l V m 



l^| < 



7T 



< 



< 



< 



2n x 


/x 4 + T N+1 








/a; 4 + t n+1 







2he~ TN + 1 +2h ^ e" 

771=^ + 2 



27rWa; 4 + r 



2fe~ T *+ 1 + 2 



df 



2/ie r «+ 1 



N+l 



Q ' N + l 



{2hr N+x + p _ T 2 +i 



27rrAr + i A /a: 4 + t 4 



Af+1 



To arrive at the last line we have used that, for x > 

-t 2 



e _t dt 



x 



e " , c 
— dt < - 



(45) 



The choice of h we make is designed to approximately equalise Ah (x) and 
this bound on Tjy. We choose h so that H = n/h = tm+i = (N + l/2)h, 
i.e., we make the choice h = sJir/(N + 1/2) given by (f2l)]). in which case 
tn+i = An = \J(N + l/2)7r, and = tfe, where is defined by ([7]). With 
this choice of h it holds that 



and that 



E N {x) = F(x) - F N (x) = e* h + T N 

\ Tn \<— + 

2nA N ^/x i + A% 
Thus we arrive at our main pointwise error bound, that 

(2tt + 1)|x| 



\E N (x)\ < m (x) := A h (\x\) 



2nA NV / x i + Al 



(46) 



N 



with h = yJ-Kj (N + 1/2) so that H = ir/h = An- We have shown this bound 
for x > 0, but the symmetries ((16)) and (fl~7|) imply that En{—x) = —En{x), 
so that (|46|) holds also for x < 0, and, by continuity, also for x = (and in 
fact En(0) = ?7Ar(0) = 0). Explicitly, for this choice of h we have that 



— A 2 

xc A N 



^{A%-x*/2) (1-e- 2 ^) 



< 



V2 



A h {x) = < 



Axe N 



(l + 20F< 



y/irA N (A N +x/y/2) (l 



-2,i-; 



\A N < 



71 



: (x 2 /2-A 2 n )(1-c- 2A n) l 



V2 



(47) 
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We will compare |£ , g(x)| to the upper bound 779(2;) in Figure [5] below. 

Note the factor cxp(— A 2 N ) = cxp(— tt(N + 1/2)) in each of the terms on 
the right hand side of (|4"T)) . so that implies that Fn(x) is exponentially 
convergent as N — > 00 for each x. Note also that Ah(x) is increasing on 
[0, jV^An) and decreasing on [|v2 An, 00), behaving asymptotically like 
kjyx^ 1 as x — > 00, where = 2e _ " 4jv /(y / 7r (1 — c~ 2j4jv )). Further, where 
^^(Iv^ A N ) denotes the limiting value of Ah(x) as x — > f \/2 from be- 
low, since 2Aj r 1 > e~ A ™, 

9y/TT A N (1 — e 2j4 «) 1 — e 5j4 iv/ 2 

Similarly, x/\^(a;) is increasing on [0, jV2An) and decreasing on [jy/2 An, 00), 
approaching the limit fcj\r as a; — > 00. Thus 

A h {\x\) < A h (^V2A N ) and xA h (\x\) < \y/2A N A h (^V2A N ) , for xeR. 

(48) 

Moreover, 



; 1 1 < -= and ; < 1, for x e K. (49) 

V^ 4 + A% ~ s/2A N v 7 ^ 4 + A % 

Combining (gHJ), ([38]). and ([4"9]>. we see that 



|F(x) - Fjv(x)| = |£?jv(a:)| < w(x) < c N , for x £ R, (50) 

ViV+1/2 

where 

(2tt + 1) e" A ? 



c N =e nN y/N+\/2 



A h {\V2A N ) 



2V2tt ^4 2 



iV 



is given by (fT5)l . Note that cjv decreases as iV increases, with c\ ~ 0.8249 
and lmi/v^ooCAT ~ 0.2080 given by (fT4l) . Of course, as observed in the in- 
troduction, this simple, explicit bound shows exponential convergence of the 
approximation Fn to F, uniformly on the real line. 

In the appendix it is shown that |F(x)| > (2 + 2^/ttx)^ 1 for x > 0, and 
that \F(x)\ > 1/2 for x < 0. Combining these bounds with (gBJ), (gSJ), and 
dHI), we see that 

\F(x)\ ~ \F(x)\ ~ \ 2c N - w ==, for x < 0, 
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where 



c* N = 2e* N 



(2tt + 1) e" A « / 1 



2tt A n \V2A n 
10V2(A + bV2^A N ) (l + 20Fe-' 3A «) ( 2 7r + l) / 1 



9y/^e*/ 2 A N (l-e- 2A2 N) ne*/ 2 A N \y/2A 



Note that c* N decreases as N increases, with c* ~ 10.4 and limjv_ i . 00 c* N = 
100e _7r / 2 /9 m 2.3. The bound (f5"Tj) shows exponential convergence of the rel- 
ative error, |i*jv(x) — F(x)\/\F(x)\, uniformly on the real line, in particular 
showing that the relative error is < 1.6 x 10~ 16 on the whole real line if 
N = 12 (see Figure [1] below). 

The above estimates use ([4"6"| and (j4"T|) to bound the maximum absolute and 
relative errors in the approximation Fn(x). Let us note that these inequalities, 
additionally, imply that Fn(x) is particularly accurate for \x\ small. For \x\ < 
A N /V2 = ^(N + 1/2)tt/2, it follows from (g6j) and (gZj that 



where 



\F(x) - F N (x)\ < r,(x) < cnM^—^ (52) 



- = ? (27T+1) 

37r 3/2 e V2 (1 - e- 2A «) w*e"/ 2 A N ' 1 J 



Note that cjv decreases as N increases, with. c\ ~ 0.17 smcl lim^v— >co c^v — 
8/(3tt 3 /V/ 2 ) w 0.10. 

In SJT] we have made claims regarding the analyticity of the approximation 
Fpj(x), considered as a function of x in the complex plane. We justify these 
claims now. One attractive feature of the modified trapezium rule approxima- 
tion I£ is that, in contrast to Ih, it is entire as a function of x. This is not 
immediately obvious: It = Ih + PCh, and PCh has simple pole singularities at 
x = e -I7r / 4 Tfc, fceZ. But Ih also has simple poles at the same points and it is 
an easy calculation to see that the residues add to zero, so that the singulari- 
ties cancel out. Since Fm{x) = It — TV, with h given by (|2l)|) . it follows that 
the singularities of Fjq{x) are those of Tjv, i.e., simple poles at ±e~ I7r / 4 £fc, for 
k = N+l, N + 2, .... Thus F^(x) is a meromorphic function and, in particular, 
is analytic in the strip |Im(x)| < Ajy/^/2 and in the first and third quadrants 
of the complex plane. 

We will note two consequences of this analyticity and the bounds that 
we have already proved. In these arguments we will use an extension of the 
maximum principle for analytic functions to unbounded domains, that if w(z) 
is analytic in an open quadrant in the complex plane, let us say Q = {z 6 C : 
< | arg(z)| < 7r/2}, and is continuous and bounded in its closure, then 



sup |u;(z)| < sup |iu(z)|, 

zeQ z£dQ 



(54) 
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where dQ denotes the boundary of the quadrant. (This sort of extension of the 
maximum principle to unbounded domains is due to Phragmen and Lindclof; 
see, e.g., [25].) 

The first consequence is that, from (fTTjl . (JT2J) , and (fT8|) . it follows that 
(TT2"]) holds if x is real or if x = \Y with Fgl, i.e., the bound p2[) holds on 
both the real and imaginary axes. Further, from (|2T))) and the asymptotics of 
erfc(x) in the complex plane [5J (7.1.23)], it follows that F(x) — > 0, uniformly 
in arg(x), for < arg(x) < 7r/2; moreover, it is clear from ((5J) that the same 
holds for .Fjv(x) and hence for En(x). Thus ([54| implies that ([12"]) holds for 
< arg(x) < 7r/2, and and (TTTJ) then imply that (JT^J) holds also for 
7r < arg(x) < 37r/4. 

It is clear from the derivations above that, if h is given by (J2SJ), then It 
also satisfies the bound (1121). i.e., 



|-F(g) - It\ < c N - (55) 
^/iV + 1/2 

this holding in the first instance for real x, then for imaginary x, and finally for 
all x in the first and third quadrants. The bound (|12j) cannot hold in the second 
or fourth quadrant because En(x) = F(x) — Fn(x) has poles there. This issue 
does not hold for F(x) — It, which is an entire function, but (|55[) cannot hold 
in the whole complex plane because this, by Liouville's theorem (|25j). would 
imply that F(x) — 1^ is a constant. What docs hold is that c~ lx (F(x) — /,*) 
is bounded in the second and fourth quadrants, this a consequence of the 
definition of 1^ and the asymptotics of e z erfc(z) at infinity. Thus it follows 
from (|54[) , and since \e~ lx \ = 1 if x is real or pure imaginary, that 

e" Tjv 

\F(x) -I* h \< c N c- XY , (56) 

y iV + 1/2 

for x = X + iY in the second and fourth quadrants. 

We can use the bound ([56| to obtain a bound on i?/v (x) in the second and 
fourth quadrants. Clearly, where Tjy is defined by (|44)l . with /i given by (j26j) . 
for a; = X + iK in the second and fourth quadrants, 

-77 N 

\F(x) - F N (x)\ < c N c~ XY +\T N \. 

yiv +1/2 

Further, arguing as below (O, if 1*1 < ^4^/(2^2) so that 

which implies that \x 2 + it 2 k \ > 1x1^/(2^2), then 

(27r + l)V2 A , _ y^r + l) 

1 W| - tt^ 7r3/ 2 exp(7r/2)(^ + l/2) ' 
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Thus, for x = X+iY in the second and fourth quadrants with \Y\ < An / (2\/2) , 



-ivN 

\F(x) - F N (x)\ < c N e~ XY (57) 

yj iV + Yj I 



where 



c N :=c N -\ - ; ; (58 

7r 3 / 2 cxp(7r/2)^7V + l/2 

The sequence cat is decreasing with c\ ~ 1.14 and Iimjv->oo cjq — limjv^oo cn ~ 
0.208. 

We observe above that the bound (|T2|) on En(x) = F(x) — Fn(x) holds for 
all complex x in the first and third quadrants of the complex plane, and on the 
boundaries of those quadrants, the real and imaginary axes, while the bound 
([57)1 holds in the second and fourth quadrants for |Im(x)| < Ajv/(2v / 2). A 
significant implication of these bounds is that they imply that the coefficients 
in the Maclaurin series of Fn(x) are close to those of F(x). Precisely, at least 
for Id < A N /V2, 



OO 

n=0 n=0 



F(x) = ^ a n x n and F N (x) = b n 



witha„ = fW(0)/n!,i) n = F { " l) {0)/n\. Thus, where M N = sup,, /— \E N (z)\, 

it follows from Cauchy's estimate J25) Theorem 10.26] and the bounds (|T2l) and 
(|B7|) that, for N > 4 so that Ajv/(2\/2) > y/n/2, 

K _ M . !*£M < M „ (if 2 < «„ (iV 2 (59) 



,V ^ + 1/2 
We will use this bound to derive (|69[) and ([70)) below 



3 Approximating C(x) and 5(a;) 

From ^ we see that, for x real, 

C(a:) = Re (\/2c iT/4 (i - F{^j2x))\ , S{x) = Im (\/2c i7r / 4 (± - F{^V/2x)) 

(60) 

Clearly, given the approximation Fn(x) to F(x), these relationships can be 
used to generate approximations for the Fresnels integrals C(x) and S(x). 
These approximations are defined, for x S K, by 

C N {x) = Rc(v / 2c-/ 4 (i-F w (VV22;))) , 

> < (61) 

S N (x) = Im (V2e™/4(1 _ F N (^/2x))) , 
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and are given explicitly in © and © ■ We note the similarity between © and 
© and the formulae Q] (7.5.3)-(7.5.4)] 

C(x) = i + f(x) sin (^ttx 2 ) - g(x) cos (^ttx 2 ) , (62) 
S{x) = | - /(x) cos (^ttx 2 ) - ,g(x) sin Qttx 2 ) , (63) 

which express C(x) and S(x) in terms of the auxiliary functions, f(x) and 
g(x), for the Fresnel integrals [TJ §7.2(iv)]. Indeed, it follows from [TJ (7.7.10)- 
(7.7.11)] that, for x > 0, /(x) and g(x) have the integral representations 

tl \ 0Fx 3 f°° e"* 2 , , . , x p i 2 e"* 2 

tlx) = — - — / o dt and man = — = / t, dt. 

2 Jo (fx 2 ) 2 +£4 (fx 2 ) 2 +t 4 

Recalling that A^r is linked to the quadrature step-size through (|2!)|) . it is clear 
that, for x > 0, y/n xa^ (fx 2 ) /Aiv an d ^/w xbjy (§^ 2 ) /^4jv can be viewed as 
quadrature approximations to these integrals. 

Table shows the Matlab code implementing ([5J and © that we use for 
the computations in the next section. Some comments of explanation are in 
order regarding the evaluation of 

sinhfisint . . 

; , (64 

cosh t + cos t 

with t = y/wAj^x, in ([8]) and ((9]). An issue in floating point arithmetic for 
evaluation of (|?54")) for larger values of \t\ is overflow; for t > 711 both sinht 
and cosht evaluate in the IEEE double precision arithmetic implemented in 
Matlab as Inf , the IEEE representation for +oo, so that (|M|) is undefined. 
But, since u is indistinguishable from u ± v in IEEE double precision real 
arithmetic if \v\/\u\ < e/2, where e = 2~ 52 w 2.22 x 10" 16 (eps in Matlab) 
is the smallest number such that 1 + e is distinguishable from 1, it makes 
sense to replace (|6"4"]) by ±1 well before this point. Precisely, for \t\ > 39, 
the expressions cosh(t) + cos(t) and exp(£)/2 have the same value in double 
precision arithmetic, as do the expressions sinhtzbsint and sign(t) exp(t)/2. 
Thus (|64|) evaluates as sign(t) in double precision arithmetic for 39 < \t\ ^ 710. 
It makes sense then, both to avoid overflow and to reduce computation time, to 
evaluate (ffi"4")) as sign(t) for \t\ > 39 (which corresponds to |x| > 39/(v / 7t^4tv) ~ 
3.51 for N = 12). 

For small t there is an additional issue in evaluation of (|64|) when the 
negative sign is chosen, that sinht ps sint rj t for small t, so that there is loss 
of precision in evaluating sinht — sint for |t| small. This is avoided in the code 
in Table by using the power series sinht — sint = 2t 3 /3! + 2t 7 /7! + . . . for 
|t| < 1, truncating this after four terms as the fifth term is negligible in double 
precision arithmetic for |t| < 1. 

The approximations © and © inherit the accuracy of Fjy(x) on the real 
line: from (|6"0"|) and (ICTI) we see that the bounds (fT5")) hold for x G K, where 
En{x) = F(x) — Fn(x). Thus the error bounds of the previous section can be 
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function [C,S] = f resnelCS(x.N) 

"/, Evaluates approximations to the Fresnel integrals C(x) and S(x). 
% x is a real scalar or matrix, 

'/, N is a positive integer controlling accuracy (suggest N=12) , 

'/, C and S are the scalars/matrices of the same size as x approximating C(x) and S(x). 
h = sqrt(pi/(N+0.5)) ; 

t = h*((N:-l:l)-0.5); AN = pi/h; rootpi = sqrt(pi); 
t2 = t.*t; t4 = t2.*t2; et2 = exp(-t2); 
x2pi2 = (pi/2)*x.*x; x4 = x2pi2 . *x2pi2; 
a = et2(l) ./(x4+t4(l)) ; b = t2(l)*a; 
for n = 2:N 

term = et2(n) ./(x4+t4(n)) ; 

a = a + term; b = b + t2(n)*term; 

end 

a = a.*x2pi2; 

mx = (rootpi*AN)*x; Mx = (rootpi/AN)*x; 
Chalf = 0.5*sign(mx) ; Shalf = Chalf ; 
select = abs(mx)<39; 
if any(select) 

mxs = mx(select); shx = sinh(mxs); sx = sin(mxs) ; 

den = . 5 . / (cos(mxs)+cosh(mxs) ) ; 

Chalf (select) = (shx+sx) . *den; 

ssdiff = shx-sx; 

select2 = abs(mxs)<l; 

if any(select2) 

mxs = mxs (select2) ; mxs3 = mxs . *mxs . *mxs ; mxs4 = mxs3.*mxs; 
ssdiff (select2) =mxs3.*(l/3 + mxs4.*( 1/2520 ... 

+ mxs4.*((l/19958400)+(0.001/653837184)*mxs4))) ; 

end 

Shalf (select) = ssdiff. *den; 

end 

cx2 = cos(x2pi2); sx2 = sin(x2pi2) ; 

C = Chalf + Mx.*(a.*sx2-b.*cx2) ; S = Shalf - Mx. * (a. *cx2+b . *sx2) ; 



Table 2 Matlab code to evaluate Cjv(a;) and Sn{%) given by (0 and (5). See ij3]for details. 



applied. In particular, from (|5U)) and (|5"!Zf it follows that both \C(x) — Cjq{x)\ 
and \S(x) — Sjq(x)\ are 

< 2c N ; for x e R, (65) 

and 

<V^c N \x\^—j, for \x\ < y/N + 1/2. (66) 

Here cn < 0.83 and cjv < 0.18 are the decreasing sequences of positive numbers 
defined by flSJ and (|55|) . respectively. 

These bounds show that Cm{x) and Sn{x) are exponentially convergent as 
N — >• oo, uniformly on the real line, so that very accurate approximations can 
be obtained with very small values of N ( ([SB]) shows that both \Cn(x) — C(x)\ 
and \S N {x)-S(x)\ are < 1.4x10" 16 on the real line for N > 11). In g] we will 
confirm the effectiveness of these approximations by numerical experiments, 
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checking the accuracy of (JSJ and by comparison with the power series [TJ 
§7.6(i)] 

00 (-1)™ (^7r) 2 ™ X in+1 °° (-1)™ (^7t) 2 ™ +1 X 4 "+ 3 

C(s)= S (2n)!(4n+l) ' ^ = S (2n + l)!(4n + 3) " (67) 

It follows from the analyticity of Fn(x) in the complex plane, discussed in 
that Fn(x) has a Maclaurin series convergent in |x| < An/V2, and from 
(f6"Tj) that Cn{x) and Sn{x) have convergent Maclaurin series representations 
in |x| < An/t/it. From the observations below (TTOJ) it is clear that, echoing 
(|6"7|) , these take the form 

00 00 
Cat(.t) = ^c na; 4 " +1 , S N (x) = ]T S „x 4 "+ 3 . (68) 

n=0 n=0 

Further, it follows from and ([59")) that the coefficients c n and s n are close 
to the corresponding coefficients of C(x) and S(x), with the difference having 
absolute value 

e -7T(iV-l/2) 

for iV > 4, where < C4 < 0.77 is the decreasing sequence of positive 
numbers given by (|58|) . This implies that, near zero, where C(x) has a simple 
zero and S(x) a zero of order three, the approximations CV(x) and Sn(x) 
retain small relative error. For Cn(x) this follows already from but to see 
this for Sn{x) we need the stronger bound implied by that, for |a;| < 1, 

5(x) - S N (x)\ <V2c N , > x 411+3 = I 1 . . 

(70) 



3.1 Other Methods for Computing Frcsncl Integrals 

Naturally, there exist already a number of effective schemes for computation of 
Frcsncl integrals, and we briefly summarise now the best of these. An effective 
computational method for smaller values of |x| is to make use of the power 
series (|67[) . These converge for all x, and very rapidly for smaller x, and so 
are widely used for computation. For example, the algorithm in the standard 
reference [23] uses these power series for |x| < 1.5. For this range, after the 
first two terms, these series are alternating series of monotonically decreasing 
terms, and the error in truncation has magnitude smaller than the first ne- 
glected term. Thus, for |x < 1.5, the errors in computing C(x) and S(x) by 
these power series truncated to N terms are < 2 x 10~ 16 and < 2.3 x 10~ 17 , 
respectively, for N = 14. 

For |x| > 1.5, [Uj recommends computation using (JBDJ) and (|2"0"|) and the 
continued fraction representation for e z erfc(z) = w(iz) given as [TJ (7.9.2)]. 
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Methods for evaluation of w(z) based on continued fraction representations 
for larger complex z (which can be applied to evaluate F(x) and hence C{x) 
and S{x)) are also discussed in Gautschi [11] and are finely tuned, to form 
TOMS "Algorithm 680", in Poppe and Wijers [2T1I32]. which achieves relative 
errors of 10 -14 over "nearly all" the complex plane by using Taylor expansions 
of degree up to 20 in an ellipse around the origin, convergents of up to order 
20 of continued fractions outside a larger ellipse, and a more expensive mix of 
Taylor expansion and continued fraction calculations in between. 

Wcideman [29| presents an alternative method of computation (the deriva- 
tion starts from the integral representation (|21l) ) which approximates w(z) by 
the polynomial 



in the transformed variable Z = (L + iz)/(L— iz). Here L = y Mj\[2 and the 
coefficients a n can be viewed as Fourier coefficients and efficiently computed 
by the FFT. We will see in Sf4]that a polynomial degree M = 36 in (fTTj) suffices 
to compute F(x) = e lx w(e 17T ^ 4 x)/2 with relative error < 10~ 15 uniformly on 
the real line. Weideman [29] argues carefully and persuasively that, in terms of 
operation counts, the work required to compute w{z) with the 10 -14 relative 
accuracy of Algorithm 680 [32] is much smaller using the approximation ([7T]) 
for intermediate values of \z\ (values in approximately the range 1.5 < \z\ < 5 
for the case arg(z) = 7r/4 which we require). 

All these approximations described above are polynomial or rational ap- 
proximations (or piecewise polynomial/rational approximations, proposing dif- 
ferent approximations on different regions). Many other authors describe ap- 
proximations of these types for computing the Fresnel integrals specifically 
with real arguments. The best of these in terms of accuracy is Cody [5] , where 
numerical coefficient values are given for piecewise rational approximations to 
C{x) and S(x) for < x < 1.6, and for piecewise rational approximations to 
f(x) and g(x) in (|62[) and (|63[) . for x > 1.6. These approximations, in their 
respective regions of validity, achieve relative errors < l() -15 - 58 pa 2.7 x 10~ 16 , 
this using rational approximations which are ratios of polynomials of degree 
< 6; in total five different approximations are used on different subintervals of 
the positive real axis. Single rational approximations, based on a "polar" ver- 
sion of ([62]) and (|63|) , are computed in [14] , but these are of limited accuracy 
(absolute errors < 4 x 10~ 8 ). 

4 Numerical Results and Comparison of Methods 

In this section we show the results of numerical computations that confirm 
and illustrate the theoretical error bounds in ^5] and SJ3] and that explore the 
accuracy and efficiency of our new methods, through qualitative and quantita- 
tive comparisons with certain of the other computational methods described 




(71) 




in 03] 
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Fig. 1 Left hand side: maximum error, max x >o \F(x) — Fn(x)\, and its upper bound (1126 

(— ), plotted against N, where F(x) is approximated by Fw(x) := e lx W3g(e 17T ^ 4 x)/2 with 

W3g(z) given by 1171 \ (— • — •), and where F(x) is approximated by F2o(x) ( ). Right hand 

side: maximum relative error, max x >o \(F(x) — Fpf(x))/F(%)\, and its upper bound (1 5 1 It 
(— ), plotted against N, where F(x) is approximated as on the left hand side. 



Turning first to the approximation (JS|) to F(x), in Figure [T] we plot against 
N the maximum values of the absolute and relative errors, \En(x)\ = \F(x) — 
Fn(x)\ and \F(x) — Fjv(x)|/|J r (a;)|, on x > 0, approximating these maximum 
values on [0,oo) by computing at 40,000 equally spaced points between and 
1,000 and replacing F{x) either by F2o(x) or by Fw{x) := e lx W3G (e 17r / 4 x)/2 
with W3q(z) given by (|7Tj) . (Wc compute Fn(x) in Matlab using the code 
in Tabic [TJ and Fuu(x) by exp(i*x.~2) . *cef (exp(i*pi/4) *x,36) /2, where 
cef .m is the function in Table 1 of [35] which evaluates ([TTj) . The choice 
M = 36 in (|7ip is made because this is the smallest value which appears to 
achieve relative errors < 10~ 15 uniformly on the positive axis, and increasing 
M further in the range M < 40 appears to lead to no increase in accuracy by 
comparison with the same approximation with M = 50; see Figure [5] below.) 
We show in the same plots the upper bounds which are the right hand sides 
of (fT2|) and (|51[) . It can be seen that the exponential convergence predicted 
by the bounds (fl2|) and (jSTj) is achieved, indeed these bounds overestimate 
their respective maximum errors by at most a factor of 10. Further, with N as 
small as 12 it appears that we achieve maximum absolute and relative errors in 
Fjv(x) which are < 2.9 x 10 -16 and < 9.3 x 10 -16 , respectively; these values are 
upper bounds whichever of the two methods for approximating F(x) accurately 
is used. (We should add a note of caution here: the different approximations 
agree to high accuracy, but the accuracy of each approximation is limited, for 
large x, by the accuracy with which e lx is computed.) 

The plots in Figure [TJ in addition to shedding light on the accuracy of 
-F/v {x) , provide independent verification of the high accuracy of the approxi- 
mation ([7T|) for w(z) proposed in [55], at least for arg(z) = 7r/4 and provided M 
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M M 

Fig. 2 Left hand side: maximum error, max x >o l-F(x) — Fw(x)\, where Fw(x) := 

c lx t0jv/(c" r/ ' 4 a;)/2 with wm{z) given by II71II . plotted against M. Right hand side: same, 
but maximum relative error, max I >o \(F(x) — Fw(x)) / F(x)\, is plotted against M. In each 
plot the two curves correspond to different methods for approximating the exact value of 
F(x), either F(x) « F 20 (x) given by JSJ (-), or F(x) « Fw(x) with M = 50 (--). 



is large enough in (|71[) . Exploring this in more detail, in Figure[5]the maximum 
absolute and relative errors in the approximation Fw(x) = e lx wm (e 17r / 4 x)/2 
for F(x), with wm{z) given by (|TTj) . are plotted against M. (The maxima, as 
in FigurefTJ are taken over 40,000 equally spaced points between and 1,000.) 
In each of the plots in Figure [5] the trend is one of exponential convergence, 
but the convergence is not monotonic and is slower than that in Figure [TJ 

In Figure [3] we plot against x the absolute and relative errors in F/v (x) 
for N = 9. On the same graphs we plot the upper bounds t]n(x) and 2(1 + 
\pR x)vn(x), respectively, with t]n(x) defined by (|4"o| . We see that the the- 
oretical error bounds are upper bounds as claimed, and that these bounds 
appear to capture the ^-dependence of the errors fairly well, for example that 
En{x) = 0{x) as x -> 0, = 0(x~ 1 ) as x — > oo, and that Em{x) reaches a 
maximum at about x = \[2A^ = y/ tt(2N + 1) (w 7.7 when N = 9). 

The above figures explore the accuracy of the approximation Fn(x). Let 
us comment now on efficiency. Most straightforward is a comparison of the 
Matlab function F(x,N) in Table [T] with computation of F(x) via the Mat- 
lab code exp(i*x . ~2) . *cef (exp(i*pi/4) *x ,36)/2 that uses cef .m from [2"9"] 
which implements (jTTj) . Both F(x,N) and cef (x,M) are optimised for efficiency 
when x is a large vector. Assuming that the time for computation in cef of the 
coefficients a n in ([TTj) is negligible, the main cost in computation of F(x) via 
cef when a; is a large vector is a complex vector exponential (for c lx ), slightly 
more than M complex vector multiplications and M additions, and 2 complex 
vector divisions (all vector operations componentwise). The major part of this 
computation is that required to evaluate the polynomial (|71|) of degree M us- 
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Fig. 3 Left hand side: absolute error, \F(x) — Fn(x)\ (— ), and its upper bound t)n(x) given 

by H46H ( ), plotted against x. Right hand side: relative error, \F(x) — Fpf(x)\/\F(x)\ (—), 

and its upper bound 2(1 + y/n x)r]N (x) ( ), plotted against x. In both plots N = 9 and 

F(x) is approximated by F2o(x). 



ing Horner's algorithm. In comparison, evaluation of F(x) using the function 
F(x,N) in Table [T] requires 2 complex vector exponentials, one complex vec- 
tor division, and slightly more than TV real vector multiplications/divisions, 
real vector additions, complex vector multiplications, and complex vector ad- 
ditions. From Figures [T]and [2] we read off that to achieve absolute and relative 
errors below 10~ 8 requires TV = 6 and M = 18; to achieve errors below 10 -15 
requires TV = 12 and M = 36. Thus it seems clear that computing F(x) via 
F(x,N) requires a substantially lower operation count than computing via cef. 
(We note, moreover, as discussed in H3.ll and in §7 of that, at least for 
intermediate values of x (1.5 < x < 5), the operation counts via cef are lower 
than those required via the method for w(z) of [2T1I22] .) 

To test whether F(x,N) is faster we have compared computation times in 
Matlab (version 7.8.0.347 (R2009a), running on a laptop with dual 2.4GHz 
P8600 Intel processors) between exp(i*x . ~2) . *cef (exp(i*pi/4) *x,36) /2 
and F(x, 12) when x is a length 10 7 vector of equally spaced numbers between 
and 1,000. The elapsed times (average of 10 executions) were 11.1 and 15.6 
seconds, respectively, so that F(x, 12) is a little less than 50% faster. 

Turning to C(x) and S(x), these can of course be computed using F(x,N) 
to calculate Fjv(i), and then using (|6ip which incurs negligible additional com- 
putation. This is entirely satisfactory except for small x, where this method 
fails to maintain small relative errors. As discussed in ^3l the Matlab func- 
tion f resnelCS .m in Table [2] directly implements ([SJ and ([9]), taking care 
in the evaluation of sinh — sin in (jTJ]) so as to achieve the high accuracy of 
Sn(x) for small |x| predicted in ([70)1 . To test the efficiency and accuracy of 
the implementation in Table [5] we have compared evaluation of C\i (x) and 
S\i{x) via f resnelCS with their evaluation via F(x,12) and (|61|) . computing 
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Ci2(x) and Si 2 (a;) at 10 7 equally spaced ^-values between and 20. The values 
of C\2(x) and Si2(x) computed by these slightly different methods differ by 
< 4.5 x 10~ 15 ; this good but not perfect agreement is because there is a dif- 
ference between exp(i(y / 7r/2 x) 2 ) and exp(i7rx/2) in floating point arithmetic. 
In this test f resnelCS requires only 67% of the computation time of comput- 
ing via F(x,12), this because the real arithmetic in f resnelCS is faster and 
because the expressions with t = \/2Ajyx, are evaluated (efficiently and 
accurately) in f resnelCS as sign(f) when |t| > 39 (corresponding to x > 3.51 
for N = 12), as discussed in 




Fig. 4 Left hand side: maximum values of \Cn(x) — C(x)| and \Sm{x) — S(x)\ on < x < 20. 
Right hand side: maximum values of \Cn(x) — C(x)\/C(x) and \Sn(x) — S(x)\/S(x) on 
< x < 20. 



These small absolute errors in Cn{x) and Sn(x), evaluated by f resnelCS, 
do not guarantee small relative errors near the only zero of C(x) and S(x) 
at x = 0. Near zero, from (|rJ71) . C(x) sa x and S(x) ~ irx 3 /6, so very ac- 
curate calculations are needed to maintain small relative errors. The bounds 
(|66[) and (|70p do, in fact, guarantee small relative errors near zero in infi- 
nite precision arithmetic, moreover predicting that these relative errors should 
decrease approximately in proportion to e -71 ^ as N increases. To see if this 
convergence is achieved in floating point arithmetic by f resnelCS, in Fig- 
ure 2] we have plotted the maximum values, on (0, 20], of the absolute errors, 
\Cn(x)—C(x)\ and \Sn(x)—S(x)\, and the relative errors, \Cn(x)—C(x)\/C(x) 
and \Sn(x) — S(x)\/S(x). In this figure we compute Sn(x) and Cjv(x) by 
fresnelCS and approximate C(x) and S(x) by 6*20(2;) and S2o(x) for x > 1.5. 
For < x < 1.5 (as recommended in [33]) we approximate by the power se- 
ries (|67[) truncated after 15 terms, and evaluating the resulting polynomials 
by the usual Horner algorithm. Exponential convergence in both the absolute 
and relative errors is seen in Figure El but we note that, while the absolute 
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errors arc < 4.5 x 10 16 for N > 11, the maximum relative error in Cn{x) is 
« 3.6 x 1CT 15 for TV = 11 and that in Sjv(a;) as large as 2.7 x 1CT 13 . These 
errors may be entirely acceptable, but the truncated power series (|67| must 
achieve smaller errors for small x, and may be cheaper to evaluate. (In fact, 
evaluating at 10 7 equally spaced points between and 1.5 takes 2.9 times 
longer in Matlab with f resnelCS than evaluating 15 terms of both the series 
(|6"T|) via Horner's algorithm.) 



5 Extensions and Concluding Remarks 

To conclude, we have presented in this paper new approximations for the 
Fresnel integrals, derived from and inspired by modified trapesium rule ap- 
proximations previously suggested for the complementary error function of 
complex argument in |18U15j . These approximations are simple to implement 
(Matlab codes are included in Tables [1] and [2J : the computation of Fn(x) re- 
quires a couple of complex exponentiations and a short summation to compute 
a quadrature sum, and that of Cjv(x) and Sn(x) evaluation of trigonomet- 
ric and hyperbolic functions and a similar short summation. The numerical 
methods are proven to converge exponentially (in absolute and relative error) , 
approximately in proportion to cxp(— ttN) where N is the number of quadra- 
ture points used. Simple explicit error bounds are provided, and the predicted 
exponential convergence is precisely observed in practice. The approximation 
Fn(x) with N = 12 quadrature points achieves close to double precision ma- 
chine precision uniformly on the real line (with the proviso that this precision 
is necessarily limited by the accuracy with which e lx can be calculated). The 
approximations Cjv(x) and Sn(x) with N = 11 are similarly accurate, except 
that the relative error in Sn(x) increases to 2.7 x 10~ 13 near x = where S(x) 
has a zero of order 3. 

Operation counts and timings carried out suggest that Fn(x) with N = 12 
may be faster than previous methods, at least for intermediate values of |a;|. 
In particular, the Matlab function in Table [1] outperforms that in Table 1 of 
[2"9"] for this application. The code for Sn(x) and Cjy(x) is faster still, but the 
power series (|67[) . truncated after 15 terms, are more accurate and efficient on 
the interval [0, 1.5], this conclusion endorsing recommendations in [23] . 

Part of the motivation for this paper was a remark in Weideman [35] re- 
garding the modified trapezium rule methods of fTMT5"] for computing erfc(z), 
that they are "very accurate, provided for given z and N [the finite number of 
quadrature points retained] the optimal stepsize h is selected. It is not easy, 
however, to determine this optimal h a priori." At least as far as computing 
erfc(z) for arg(z) = — 7r/4 is concerned (which, by ([20]) . is the same as com- 
puting F(x)) this problem is solved in this paper, so that the effectiveness 
of the modified trapezium rule methods of [TBI|15[|2"9"] is clearly demonstrated. 
We hope that the methodology and positive results of this paper will inspire 
further applications of this truncated, modified trapezium rule method. 
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With respect to this hope, most obviously the results in this paper suggest 
a revisit of the methods of [TMT5) for erfc(z). Clearly, ([2Tj)l suggests erfcjv(z) := 
2Fjs[(e 17T / 4 z), given explicitly as 

crfcAr(z) = 7 + -^c > — — (72) 

V > c 2A N z + l £ z 2 + t 2 



k=l 



A- 



as an approximation for erfc(z). (For < Re(z) < ^4 at, this is precisely the 
approximation of |15) truncated to N quadrature points and with the par- 
ticular choice (|26p for h made.) The results of show that (fT2"]) holds for 
< arg(x) < 7r/2 and for 7r < arg(s) < 37r/4 which implies that 

-ttN -ttN 

\eric(z) - erfcjv(z)| < 2c N ; < 2 - (73) 

y/N + 1/2 v/iV+1/2' 

for |arg(z)| < 7r/4 and 37r/4 < arg(z) < 57r/4. This is a strong result in 
37r/4 < arg(z) < 57r/4, where it is known that |erfc(z)| > 1 so that ([73^) is 
a bound on both the absolute and relative error. However, in | arg(z)| < 7r/4 
the bound ([75)) is less satisfactory. In particular, since erfc(x) ~ c~ x /(y/wx) 
as x — > +00, for larger x > (J73J) docs not guarantee small relative errors. 
Indeed, erfcjv^) ~ 2c~ 2AnX has the wrong asymptotic behaviour as x — >• +oo. 

A large part of a possible fix and analysis is already in [T5] and [TS] (see 
equations (7)- (8) in [T5] and cf. (f2~T|) . [T5]l namely to discard the first term in 
([73)) for Re(z) > Ajy, so that erfc(z) is approximated by 

2z 2 N e~*fc 
ertc' N (z) = R N (z) + -i- e"* 2 ^ -J— j , (74) 



fe=l 



where 

R N (z) 



2/(e 2ANZ + 1), Rc(z) < A 



N- 



0, Rc(z) > A N . 

Figure [5] plots the supremum, on < x < 25, of the absolute and relative errors 
in erfc^r(x) against N, computing erfc(x) with the inbuilt Matlab function 
erfc, and with erfc^o- Clearly, both plots show exponential convergence at 
a rate approximately proportional to e~ vN . The absolute error in erfc' w is 
< 4.5 x 10~ 16 for N > 10 and the relative error < 6.7 x 1CT 16 for N = 12, 
while the maximum relative error of the standard Matlab function is limited 
to about 5.7 x 10 -14 . We have not computed any theoretical upper bounds for 
these errors, but the fact that erfc^x) is discontinuous at An implies lower 
bounds: in particular, the supremum of the relative error in any neighbourhood 
of An must be > l/(erfc(AAr)(exp(2A^) + 1)), i.e., half the discontinuity at 
An divided by erfc(-Ajv). From lf4"5)) , erfc(ar) < e~ x /(y/irx) for x > 0, so we 
see that 

gup |crf4(x)-crfc (a; )| ^ ^ A N c^ ^ ^^-Aj (?5) 
o<x<i+a n erfc(x) e 2A N + 1 
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In the right hand figure we plot this lower bound which accurately predicts the 
maximum error, this suggestive that the small size of error present is associated 
with the discontinuity in erfc^. 




4 6 8 10 12 4 6 8 10 12 

N N 

Fig. 5 Left hand side: supremum of \erfc' N (x) — erfc(x)| on < x < 25 and e — a n /100 (— .). 
Right hand side: supremum of |erfc^y(:f;) — erfc(x)|/erfc(a;) on < x < 25 and its theoretical 

lower bound, ^/tt Aj^e~ A N (— .). In each plot two of the curves differ only in how crfc(x) is 
approximated: the Matlab built in function (— ) and erfcjQ ( ). 



These are encouraging results, but further work is needed, revisiting [18U15I 
I19j . to develop a fully discrete and accurate modified trapezium rule method 
for erfc(z) for all z in the complex plane. Note, for example, that, while Figure 
[5] shows that erfc' 12 (z) has relative error close to machine precision on the 
positive real axis, on the imaginary axis this approximation is singular at 
z = ±iifc, for k > 13. 

We finish by flagging that the modified trapezium rule method that we 
have used in this paper is applicable widely to the evaluation of integrals on 
the real line of functions that are analytic but with poles near the real axis. 
Indeed, general theories of the method are presented in Bialecki [3], Hunter 
[TB] (and see [5], [Ml §5-1.4]), and in the thesis of one of the authors j!7j . 
where the emphasis is on the particular case (|23p . where the analytic function 
f(t) = 0(1) as t — > ±00. Integrals of the form ([25)1 arise in probabilistic 
applications [3] and as representations in integral form of solutions to linear 
PDEs with constant coefficients, after solution by Fourier transform methods 
and deformation of the path of integration to a steepest descent path. One 
example which continues to be the subject of computational studies p i^OllTB] 
is the Green's function for the Hclmholtz equation Au+k 2 u = in a half-space 
with an impedance boundary condition, du/dn = ikfiu. Representations for 
this Green's function in terms of a steepest descent path integral of the form 
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(T2"3"]) . in both the 2D and 3D cases, are given in [5], and the application of the 
truncated modified trapezium rule method is discussed in [17] . 
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A Appendix 



In this appendix we prove the bounds 



7^ > 1^)1 > | w * r | 

2 2 c 17r / 4 + Jn x\ 



2 VT 



> 



+ v 27r x + nx 2 



2y/TTx' 



(76) 



for x > 0. The lower bounds in (|76[) . which appear to be new, are used in $5] 
to prove an upper bound on the relative error in the approximation Fn(x) to 
F(x). From (f?6"|) and (fTB|) we immediately deduce bounds for negative argu- 
ments which are also used in that 



l>\F(-x)\>± forx>0. 



(77) 



We can also immediately deduce bounds on the version of the Fresncl integral 
defined by (j?]), and on the complementary error function of argument ±7r/4 
(via ((20)) and that erfc(z) = erfc(^)), for example that 



\T(x)\ > 



1 



V2~ 



and 



erfc 



(e^x) 



> 



1 



1 + 



for x > 0. (78) 



The remainder of this appendix is the proof of (fTB]) . But note first that 
both 

L x {x) := — 1 and L 2 (x) := - — (79) 



2Vl + V^nx + nx 2 2 

are sharp lower bounds for |.F(x)| for x = (since \F (0)| : 
i) and in the limit x — > +00 (since |-F(a;)| ~ (2y / 7ra;)~ 1 



-2Vttx 
^(0) = 



i 2 (0) = 
as x — > +00, and 

L\(x) and £2(2) have the same asymptotic behaviour). In fact, numerical 
computations (approximating |F(x)| by |Fjv(x)| with N = 12) suggest that 
1 < |F(a;)|/ii(x) < 1.25 and 1 < F(x)/L 2 (x) < 1.35, for x > 0, so that these 
are rather sharp lower bounds over the whole positive real axis. 

To prove ([76]) we note first of all that, from (f2"2")> . on substituting t = ux. 



\F{x)\ = 



1 

2^ 



e~ x u (1 - iu 2 ) 
1 + u 4 



1 

2^ 



\g (x)-i 9l (x)\ (80) 
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where, for n = 0, 1, 



9n{x) :-- 



oo _ x u „,1n 



1 +U 4 



■ du. 



Clearly, g n (x) > is well-defined for all n and all x > by this definition, 
and also for x = for n = 0, 1, with (this computation done, e.g., by contour 
integration) go(0) = 5i(0) = tt/V2. Further, for x > (and a; = for n = 0), 



so that 



9n( x ) = -2xg n +i(x) < 0, 



g n {x) < ffn(0) = — y=, for x > and n = 0, 1. 
v2 



(81) 



(82) 



Using this last inequality in ((80|) gives |-F(x)| < ^, for a; > 0. Moreover, ([80 
implies that 



c>/ 4 + 



1^(^)1 > G(x) := -^-Re ((l + + i) ( 5o (x) - i 9l (x))) 



=- ((l + \^a;) 5o (x)+ 3 i(x)) 



V2 

Clearly, ([71]) will follow if we can show that G(x) > 1 for a; > 0. 
Now, for x > 0, 



(83) 



go{x) +52(2) 



dw 



(84) 



so that, for x > 0, 



^(x) = —2xgi(x) — — 2 Vtt 4- 2xg (x) 



and 



From 



9o(x) = 5o(0) + y g' (t)dt = — -2j t 9l (t)dt, 

pX pX 

ff 1 (s) = si(o) + y o <?!(*) d* = ^=- 2^ + 2 y^ t 5o (t)<tt. 

we see that <7o(^) < tt/V2, and then from (|86l) that 
7T „ /— t 2 

Si 0*0 < ^= - 2V7T.T + -j=X . 



(85) 
(86) 



It follows from flSSI) that 



7T 7T 2 4 _ o 7T 



X , 



(87) 



and then from ([SB]) that 

. 7T ^- 7T 
5l (X) > — — 2y/TTX ~\ -= 



V2 



7= x +—VKX 

V2 2V2 15 



irV2 
12 ! 
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all these bounds holding for x > 0. Using these lower bounds in ((83")) we see 
that, for x > 0, 

V; " v/2^ 7 v^V 3/ 6 V 7T V 4 15/ 12 



1 + -7= 



where 

We will show now that ho(x) > for < x < 1 which will show that 
G(x) > 1 for < x < 1. To see this we observe that = xhi(x) where 
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X 

15 J 




7^ 




1 


4, 




3 



hAx) = -2 7T - - + — x- 2tt x 2 x 3 

n ' V 3/ 2 V 15 7 12 

< - y + 7x - Ax 2 = - \2x - 
so that h' (x) < for x > 0. Thus, for < x < 1, 

, , , , . , Sv 7 ^ 7T 2 

M*) > /*>(!) = - 2 - 15 > 0> 

We have shown that G(x) > 1 for < a; < 1. It remains to show that 
G(x) > 1 for x > 1. To see this we make use of (|83[) and ((84)) which give that, 
for x > 0, 

G(x)>-± r (l + V2^x) (^-g 2 (x) 

= 1 + ^Mi + ^)* ( 4 <89) 



A simple upper bound on g 2 (x) is 



e- x2u \ 4 du= ^J, (90) 



where to obtain this last value we integrate by parts twice and then use ((84 
Thus, for x > 0, 

V2^x(G(x) - 1) > 1 _ JL 



4x 4 4x 3 

Thus G{x) > 1 if x is large enough, in particular if x > \/2. 
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To show that G(x) > 1 for 1 < x < y/2 we need a sharper upper bound on 
g-x (x) than (|90[) . To obtain this upper bound wc write 

f°° e"* 2 t 4 
xg 2 (x)=I:= J^^-^dt, 

and approximate this integral by the trapezium rule as 

' l ^ x^ + nW ~ 2h 2^ x 4 + n 4 h 4 ■ 

n— — oo n—1 

Arguing as in §3] (or see |MJ §5.1.4]), the error in this trapezium rule approx- 
imation is / — Ih = PCh + E^, where PCh, a pole contribution, and E* h are 
given by 

PCh = 27ri(r + n) and E* h = [ f(z)(l + g(z)) dz. 

Jr H 

Here, f(z) = cT 2 ^ z 4 j {x 4 + z 4 ), g(z) = -~icot(nz/h), and, for m = 0,1, r m is 
the residue of f(z)(l + g(z)) at the pole z m , with z = xc 17T ^ 4 and z\ = iz . 
(The difference in definition of g compared to §3] is because the trapezium rule 
here is shifted by h/2 compared to §SJ note that g(z + h/2) — itan(7rz//i) as in 
§3]) Further, is the line Qz = H as in §3J and we assume that H > x/y/2, 
so that the poles at zq and z\ lie between Ph and the real axis. 
Arguing as in §3J in particular using the bound we see that 



\E* h \ < sup 



z 4 (l + g(z)) 



00 , , 9 /tt ff 4 p H —2-kH/H 



{H 4 -x 4 ){l-c-^ H / h )' 



provided that H > x. Further, 

1 ; 



'0 



£-(1 + g(zo)) = -e^ 4 - x )x(l + g(zo)), 



2(z% + ix 2 ) K ' * y v " 4 



so that \r \ = x\l + .g(z )|/4. Similarly, |n| = x\l + g(zi)|/4, so that, using 

7rTf ,-v / 2 7r2;/h 

|PC h | < ^ TOe _ , . 

1 1 - 1 _ e -V2*x/h 

Finally, using (|45p. 
Thus 

/</fc + |PG A | + |^| 



< 



c 



'^n 4 /i 4 e-^ 2 '' 2 2 7 rxe"^^/ /l 2v^ff 4 c^ 2 - 2 ^/' 1 



a; 4 -r-n 4 /i 4 AT/i l _ e -V2 -kx/H (H 4 - x 4 )(l - e -^ H / h ) ' 

n—1 v y v 7 
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provided that H > x. In particular, for 1 < x < y/2, choosing h = 1, H = 7T, 
and TV = 2, and noting that le"^ 11 is decreasing on x > 1, we see that 

2 32 1 2nxc~^ 7IX 2J^n i e-^ 
x 92{x) < i ^ + „ 4 ^4 , 1^ + — + : -7WZZ + — IZTTi 



c(x 4 + l) e 4 (a; 4 + 16) 2e 4 l _ e -V2wx (7r 4 -x 4 )(l 
where 

_32_ J_ 2^"^ 2 v We~ 7r2 
! ~ 17c 4 + 2e 4 + i _ c -v^ + - 4)(1 - c" 2 - 2 ) ~ 

Thus 

/I ^ \ /l ^ \ / 2 

+ d 



4= + V2x ) a^x) < := ( + V2~ , 

Now, for 1 < i < \/2, since x 4 /(l + cc 4 ) 2 is decreasing oni> 1, 

2 V2^ - 3V2^rx 4 - 4x 3 



(1 + x 4 )- 
2 Sv^x 4 24 



Thus, for 1 < x < V2, it follows from ([89]) that 

v/2ttx(G?(x) - 1) > 1 - ( + V2x ) x 32 (x) > 1 - 



> 1 - W(l) = 1 - f + N/2J (c^ 1 + 5) > 0. 
Thus G(x) > 1 for 1 < x < s/2 and the proof is complete. 



